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_C Abstract 

A version of the fundamental mean-square convergence theorem is proved for stochastic 
differential equations (SDE) which coefficients arc allowed to grow polynomially at infinity 
and which satisfy a one-sided Lipschitz condition. The theorem is illustrated on a number 
of particular numerical methods, including a special balanced scheme and fully implicit 
, methods. Some numerical tests are presented. 

AMS 2000 subject classification. Primary 60H35; secondary 65C30, 60H10. 
Keywords. SDEs with nonglobally Lipschitz coefficients, numerical integration of SDEs 
in the mean-square sense; balanced methods; fully implicit methods; strong convergence; 
almost sure convergence. 



X " 1 Introduction 



Let (fi,.F, P) be a probability space and {w{t)^Ff) = {{w\{t), . . . , w m (t)) T , F™) be an Tri- 
dimensional standard Wiener process, where Tf , < t < T, is an increasing family of a- 
subalgebras of T induced by w(t). We consider the system of Ito stochastic differential equations 
(SDE): 

m 

dX = a(t,X)dt + ^2a r (t,X)dw r (t), t € (t ,T], X(t ) = X , (1.1) 

r=l 

where X, a, a r are ci-dimensional column-vectors and Xq is independent of w. We suppose that 
any solution X tQj x (t) of ([Lip is regular on [io,^ 1 ]- We recall [3] that a process is called regular 
if it is defined for all to < t <T. 

In traditional numerical analysis for SDE [12|. [91 [15] it is assumed that the SDE coefficients 
are globally Lipschitz which is a significant limitation taking into account that most of the models 
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of applicable interest have coefficients which grow faster at infinity than a linear function. If 
the global Lipschitz condition is violated, the convergence of many usual numerical methods can 
disappear (see, e.g., \22\ \E[ [6| 116]). This has been the motivation for the recent interest in both 
theoretical support of existing numerical methods and developing new methods or approaches 
for solving SDE under nonglobal Lipschitz assumptions on the coefficients. 

In most of SDE applications (e.g., in molecular dynamics, financial engineering and other 
problems of mathematical physics), one is interested in simulating averages E,ip(X(T)) of the 
solution to SDE - the task for which the weak-sense SDE approximation is sufficient and effective 
[12\ [T5] . The problem with divergence of weak-sense schemes was addressed in [16] (see also [T7J ) 
for simulation of averages at finite time and also of ergodic limits when ensemble averaging is 
used. The concept of rejecting exploding trajectories proposed and justified in [16] allows us to 
use any numerical method for solving SDE with nonglobally Lipschitz coefficients for estimating 
averages. Following this concept, we do not take into account the approximate trajectories X(t) 
which leave a sufficiently large ball Sr := {x : \x\ < R} during the time T. See other approaches 
for resolving this problem in the context of computing averages, including the case of simulating 
ergodic limits via time averaging, e.g. in \22 \ I10 |, IT] . 

In this paper we deal with mean-square (strong) approximation of SDE with nonglobal Lips- 
chitz coefficients. Mean-square schemes have their own area of applicability (e.g. for simulating 
scenarios, visualization of stochastic dynamics, filtering, etc., see further discussion on this in 
[9j [T5| [7] and references therein). Furthermore, mean-square approximation is of theoretical 
interest and it also provides fundamental insight for weak-sense schemes. 

We note that in the case of weak approximation we often have to simulate large dimensional 
complicated stochastic systems using the Monte Carlo technique (or time averaging), which is 
typical for molecular dynamics applications, or we have to perform calculations on a daily basis, 
which is usual, e.g., in financial applications. Hence the cost per step of a weak numerical 
integrator should be low, which, in particular, essentially prohibits the use of implicit methods. 
In contrast, areas of applicability of mean-square schemes, as a rule, do not involve simulation of 
a large number of trajectories or over very long time periods and, consequently, there are more 
relaxed requirements on the cost per step of mean-square schemes and efficient and reliable 
implicit schemes have practical interest. There have been a number of recent works, including 
[6j [5j [HJ [2T| [7] (see also the references therein) , where strong schemes for SDE with nonglobal 
Lipschitz coefficients were considered. An extended literature review on this topic is available 
in[7J. 

In this paper we give a variant of the fundamental mean-square convergence theorem in the 

case of SDE with nonglobal Lipschitz coefficients, which is analogous to Milstein's fundamental 

theorem for the global Lipschitz case [11] (see also |12t 115]). More precisely, we assume that the 

SDE coefficients can grow polynomially at infinity and satisfy a one-sided Lipschitz condition. 

The theorem is stated in Section [2] and proved in Appendix S Its corollary on almost sure 

convergence is also given. In Section [2] we start discussion on applicability of the fundamental 
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theorem, including its application to the drift-implicit Euler scheme and thus establish its order 
of convergence. Strong convergence (but without order) of this scheme was proved for SDE 
with nonglobal Lipschitz drift and diffusion in \21\ [7] . A particular balanced method (see the 
class of balanced methods in |13| 115] ) is proposed and its convergence with order 1/2 in the 
nonglobal Lipschitz setting is proved in Section O In Section [J] we revisit fully implicit (i.e., 
implicit both in drift and diffusion) mean-square schemes proposed in p3] (see also |15|). In 
[144 [T5] their convergence was proved for SDE with globally Lipschitz coefficients. Here we relax 
these conditions as the drift is required to satisfy only a one-sided Lipschitz condition and be 
of not faster than polynomial growth at infinity. Some numerical experiments supporting our 
results are presented in Section [5j 

2 Fundamental theorem 

Let Xt 0t Xo(t) = X(t), to < t < T, be & solution of the system (jl.ip . In what follows we will 
assume the following. 

Assumption 2.1. (i) The initial condition is such that 

E|X | 2p < K < oo, for all p > 1. (2.1) 

(ii) For a sufficiently large po > 1 there exists a constant c\ > such that 

2 — 1 rn 

(x-y,a(t,x)-a(t,y))+ P ° 2 ^ \a r (t, x)-a r (t, y)\ 2 < Cl \x-y\ 2 , te[t ,T], x,yeR d . (2.2) 

r=l 

(Hi) There exist ci > and x > 1 such that 

\a(t,x)-a(t,y)\ 2 <c 2 (l + \x\ 2K - 2 + \y\ 2 *- 2 )\x-y\ 2 , t € [t ,T], x,y£R d . (2.3) 

We note that (|2.2p implies that 
2 — 3 m 

(x,a(t,x)) + -^- ^2\a r (t,x)\ 2 < co + c[\x\ 2 , t G [t , T], x G R d , (2.4) 

2 r=l 

where c = \a(t, 0)\ 2 /2 + ( 2 po-3K2 Po -i) ^ (^(^o)] 2 and c[ = c x + 1/2. The inequality (12^) 
together with (|2.ip is sufficient to ensure finiteness of moments [3]: there is K > 

E|A 40iXo (t)| 2p <A-(l + E|A | 2p ), 1 <p<po-l, *€[to,r]. (2.5) 

Also, (|2.3p implies that 

a?)| 2 < c 3 + C2|a:| 2j< , t G [to, T], x G M d , (2.6) 
where C3 = 2|a(i, 0))| 2 + 2c2(x- l)/x and c' 2 = 2c2(l + x)/x. 
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Introduce the one-step approximation X t x (t + h), to < t < t + h < T, for the solution 
Xt, x (t + ^) °f (|l-ip . which depends on the initial point (t, x), a time step h, and {wi{9) — 
wi(t), . . . , w m {6) — w m (t), t < 9 < t + h} and which is defined as follows: 

X t , x (t + h) = x + A(t,x,h-Wi(6) - Wi(t), i = l,...,m, t<6 <t + h). (2.7) 

Using the one-step approximation (j2.7|) , we recurrently construct the approximation (X k , Tt k ) , k = 
0,...,N, tk+i —tk = h k+ i , T/v = T : 

X = X(t ), X k+1 = X tkj x k (tk+i) (2.8) 
= X k +A(t k ,X k ,h k+1 ;Wi(6) - Wi(t k ), i = l,...,m, t k < 6 < t k+1 ). 

The following theorem is a generalization of Milstein's fundamental theorem (see also 
[12\ [T5"l Chapter 1]) from the global to nonglobal Lipschitz case. It also has similarities with 
a strong convergence theorem in [5] proved for the case of nonglobal Lipschitz drift, global 
Lipschitz diffusion and Euler-type schemes. 

For simplicity, we will consider a uniform time discretization, i.e. h k = h for all k. 

Theorem 2.1. Suppose 

(i) Assumption \2. 1\ holds: 

(ii) The one- step approximation Xt tX (t + h) from (|2.7p has the following orders of accuracy: 
for some p > 1 there are a > 1, ho > 0, and K > such that for arbitrary to < t < T — h, 
x G H d , and all < h < ho : 

\E[X t>x (t + h)- X t , x (t + h)]\< K(l + \x\ 2a ) l ' 2 h^ , (2.9) 

[E\X ttX (t + h)- X t>x (t + h)\ 2p ] V(2p) < K{1 + \x\ 2ap fl^h^ (2.10) 

with 

<Z2>^, qx>q2 + \\ (2.11) 

(Hi) The approximation X k from (12. 8p has finite moments, i.e., for some p > 1 there are (3 > 1, 
ho > 0, and K > suc/i i/iaf /or a// < /i < ho and all k = 0, . . . , N: 

E\X k \ 2p < K(l + E\X \ 2p P). (2.12) 

T/ien /or any N and k = 0, 1, . . . , N the following inequality holds: 

[E\X to>Xo (t k ) - X tQ , Xo (t k )\ 2p ] 1/{2p) < K(l + E\Xo\ 2 ™) 1/(2p) h q2 - 1/2 , (2.13) 

where K > and 7 > 1 do not depend on h and k, i.e., the order of accuracy of the method 
(gSD isq = q 2 - 1/2. 

The theorem is proved in Appendix 1X1 and it uses the following lemma. 
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Lemma 2.1. Suppose Assumption \2. 1\ holds. For the representation 



X t , x (t + 9)-X t} y(t + 9) 



x-y + Z t ,x, y (t + 9) 



(2.14) 



we have for 1 < p < (po — 1) jx : 



E\X t>x (t + h)- X t Jt + h)\ 2p <\x- y\ 2p (l + Kh) , 



(2.15) 



E\Z t>x>y (t + h)\ 2p <K{l + \x 



2x-2 



+ \y\ 2K - 2 ) p/2 \x-y\ 2p h p . 



(2.16) 



This lemma is proved in Appendix [B] Theorem 12.11 has the following corollary. 

Corollary 2.1. In the setting of Theorem \2. 1\ for p > l/(2q) in (j2. 13|) . there is < e < q and 
an a.s. finite random variable C(u) > such that 



i.e., the method (I2.8P for (jl.ip converges with order q — e a.s. 

The corollary is proved using the Borel-Cantelli-type of arguments (see, e.g. [2~1 ITS]). 

2.1 Discussion 

In this section we make a number of observations concerning Theorem 12.11 

1. As a rule, it is not difficult to check the conditions (|2.9p - (|2.10p following the usual routine 
calculations as in the global Lipschitz case |12t [91 [15]. We note that in order to achieve the 
optimal q\ and q2 in (|2.9p - (|2.1(jp additional assumptions on smoothness of a(t,x) and a r (t,x) 
are usually needed. 

In contrast to the conditions (|2.9p - (|2,10p . checking the condition (|2.12p on moments of a 
method Xk is often rather difficult. In the case of global Lipschitz coefficients, boundedness of 
moments of X^ is just direct implication of the boundedness of moments of the SDE solution and 
the one-step properties of the method (see |15l Lemma 1.1.5]). There is no result of this type 
in the case of nonglobal Lipschitz SDE and each scheme requires a special consideration. For 
a number of strong schemes boundedness of moments in nonglobal Lipschitz cases were proved 
(see, e.g. [6j [5j [U [22]). In Section[3]we show boundedness of moments for a balanced method 
and in Section [J] for fully implicit methods. 

Roughly speaking, Theorem 12.11 says that if moments of X^ are bounded and the scheme 
was proved to be convergent with order q in the global Lipschitz case then the scheme has the 
same convergence order q in the considered nonglobal Lipschitz case. 

2. Assumptions and the statement of Theorem 12.11 include the famous fundamental theorem 
of Milstein [11] proved under the global conditions on the SDE coefficients (of course, as discussed 
in the previous point, this case does not need the assumption (|2.12p ). Though the main focus 
here is on cases when drift and diffusion can grow faster than a linear function at infinity, we 
note that the assumptions also include the case when the diffusion coefficient grows slower than 



X to , Xo (t k ) - X k \ <C(ui)hi- £ 
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linear function at infinity, e.g. they cover so-called CIR process which is used in modelling short 
interest rates and stochastic volatility in financial engineering. 
3. Consider the drift-implicit scheme |X5|, p. 30]: 

m 

X k+l = X k + a(t k+1 ,X k+1 )h + °r(tk, X k )£ rk Vh, (2.17) 

where £ rfc = (w r (t k+ i) — w T {t k ))/Vh are Gaussian AA(0, 1) i.i.d. random variables. Assume 
that the coefficients a(t,x) and a r (t,x) have continuous first-order partial derivatives in t and 
the coefficient a(t, x) also has continuous first-order partial derivatives in x % and that all these 
derivatives and the coefficients themselves satisfy inequalities of the form ([2.3[) . It is not difficult 
to show that the one-step approximation corresponding to (|2.17p satisfies (|2.9p and (|2.10p with 
qi = 2 and qi = 1, respectively. Its boundedness of moments, in particular, under the condition 
(|2.4p for time steps h < l/(2ci), is proved in [7J. Then, due to Theorem 12. 1\ (|2.17p converges 
with mean-square order q = 1/2 (note that for q = 1/2, it is sufficient to have q± = 3/2 
which can be obtained under lesser smoothness of a). Further, in the case of additive noise (i.e., 
a r (t, x) = cr r (t), r = 1, . . . , m), q\ = 2 and qi =3/2 and (|2.17p converges with mean-square order 
1 due to Theorem 12.11 We note that convergence of (|2.17p with order 1 /2 in the global Lipschitz 
case is well known |12[ (9[ [15]; in the case of nonglobal Lipschitz drift and global Lipschitz 
diffusion was proved in [BJ [5] (see also related results in [21 [22]); and its strong convergence 
without order under Assumption 12.11 was proved in |21l [7] . 

5. Due to the bound (|2.5p on the moments of the solution X(t), it would be natural to 
require that f3 in (|2.12p should be equal to 1. Indeed, (|2.12p with (3 = 1 holds for the drift- 
implicit method (|2.17p [7J and for fully implicit methods (see Section H]). However, this is not 
the case for tamed- type methods (see [8]) or the balanced method from Section [3l 

6. The constant K in f|2. 13|) depends on p, to, T as well as on the SDE coefficients. The 
constant 7 in (|2.13p depends on a, j3 and x. 



3 A balanced method 

In this section we introduce a particular balanced scheme from the class of balanced methods 
introduced in [13] (see also [15j ) and prove its mean-square convergence with order 1/2 using 
Theorem 12. li As far as we know, this variant of balanced schemes has not been considered 
before. In Section [5] we test the balanced scheme on a model problem and demonstrate that 
it is more efficient than the tamed scheme (|5.2p (see Section [5]) from [7] . We also note that 
it was mentioned in [7] that a balanced scheme suitable for the nonglobal Lipschitz case could 
potentially be derived. 

Consider the following balanced-type scheme for (|l.ip : 

x =x a(t k ,X k )h + ^ =1 a r (t k ,X k )^ rk y/h 

k+1 k l + h\a(t k ,X k )\+y/hY^ x \*& ki X k )ti rk \' 
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where £ rk are Gaussian jV(0, 1) i.i.d. random variables. 

We prove two lemmas which show that the scheme (|3.1|) satisfies the conditions of Theo- 
rem 12.11 The first lemma is on boundedness of moments, which uses a stopping time technique 
(see also, e.g. [T6j 17]). 

Lemma 3.1. Suppose Assumption \2.1\ holds with sufficiently large pq. For all natural N and 
all k = 0, . . . , N the following inequality holds for moments of the scheme (|3,ip : 

Po-1 1 



E\X k \ 2p < K(l +E\X \ 2p P), Kp< 

~ v 1 U| ; ' ~ y ~ 4(3x-2) 2 

with some constants (3 > 1 and K > independent of h and k. 



(3.2) 



Proof. In the proof we shall use the letter K to denote various constants which are independent 
of h and k. We note in passing that the case x = 1 (i.e., when a(t,x) is globally Lipschitz) is 
trivial. 

The following elementary consequence of the inequalities (|2.4p and (j2.6|) will be used in the 
proof: for any C\ > and C2 > : 

™ ^2 

Ci^|a r (t,x)| 2 < (2C lC +-f)(l + |x| 2 ) + C7 2 |a(t,x)| 2 (3.3) 



r=l 



c 2 

<2 



< (2C lC +^)(l + |x| 2 ) + C 2C (l + |x| 2 -), 



C 

where c = max(co, c[, c' 2 , C3). 
We observe that 

l-Xfc+il < l^fcl + 1 < \X \ + (k + l). 
Let i? > be a sufficiently large number. Introduce the events 

n Rtk :={u:\Xi\<R, I = 0, . . . , k}, 

and their compliments AR tk . We first prove the lemma for integer p > 1. We have 

E X ^ fc+ »|X fc+1 | 2 P < E X{lRk (u;)\X k+1 \ 2 P = Exn^MKXfc+i - **) + X k \ 2p 
< E X n R k (u)\X k \ 2 P +E XfiRifc (w) |X fc | 2p " 2 [2p(X fc ,X fc+1 - X fc ) +p(2p- l)|X fc+1 - X fe 

Consider the second term in the right-hand side of ()3.6|) : 

EXn Hfc (w) |^ fc | 2p - 2 [2p(X fc ,X fc+1 -X fe ) +p(2p- -X fc | 2 ] 



(3.4) 



(3.5) 



(3.6) 

1 21 



(3.7) 



|2p-2 



E 



Xh 



a(t k ,X k )h + X]"Li a r (tk,X k )€ rk Vh 
1 + /i|a(i fc ,X fc )| + v^^rli Wr(tk, X k )( rk \ 



2p-l 



a(t k , X k )h + YZLi cr r(tk,X k )£ irk Vh 



l + h\a(t k ,X k )\+Vh^ =1 \a r (t k ,X k )^ rk \ 
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Since 



and for I ^ r 



E 



E 



1 + h\a(t k ,X k )\ + y/hYZL\ Wr(t k ,X k )£ rk \ 
&r{t k , X k )^ rk Vhai(t k , X k )£ lk Vh 



(i + h\a(t k ,x k )\ + VhET=i Wr{t k ,x k )i rk \y 

the conditional expectation in (|3.7p becomes 

a(t k , X k )h + Y%Li °Y(*fc, X k )£ rk Vh 



0. 



(3.8) 



(3.9) 



.4 



E 



X k , 



+ 



E 



+ 



2p-l 



i + x fc )| + VhY,?=i WAtk, x k )^ rk \ 

a(t k ,X k )h + ^r{tk,X k )^ rk \/h 



(3.10) 



1 + /l|a(tfc,X fc )| + \//lI]r=l IM^J-XfcKrfcl 

(X k ,a{t k ,X k )h) 

1 + h\a(t k ,X k )\ + \/h\ J2T=i a r{t k ,X k )£ rk \ 

2 P - 1 a 2 (t k , x k )h 2 + /i x:r=i (°r(t k ,x k )t rk ) 2 



tk 



1 + h\a(t k ,X k )\ + Z)J!LiO>(*fc,^jfc)£rfc 



< E 



(X k ,a(t k ,X k )h) 



1 + h|a(£ fc , Jffc)| + V^S™ =1 |o>(*fc,^fc)£ 



rk 



2p-l 



2e2 



^E™iM^,^)r£ 



rfc 



2 l + /i|a(4,^ fe )| + V^E^ikr(iifc,^ 



+ 



^-±a 2 (t k ,X k )h 2 



E 



+ 



(X k ,a(t k ,X k )h) + &=±hYT =1 \°r{t k ,X k )\< 
1 + h\a(t k ,X k )\ + \/^X)^Ll kr(ifc,^fc)Crfel 

2p-i hY:uwr(t k ,x k )\ 2 (e rk -i) 



2 1 + /i|a(i fc ,X fc )| + VhY17=i Wr(t k ,X k )^ rk 
Using ([23]) and (f2T6|) . we obtain 

A < Co/t + cilXfcl 2 /! 

_ 771, 

+^-^^X;k r (i fe ,x fe )i 2 E 

^ r=l 



+ 



(3.11) 



1 + h\a(t k ,X k )\ + V^X^j \a r (t k ,X k )£ 



rk 



T 



tk 
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For the expectation in the second term in (|3.1ip . we obtain 



E 



E 



1 + h\a(t k ,X k )\ + \/hYT=\ \ a r(tk,X k )£ rk 



(3.12) 



-E 



(frk ~ 1) 



1 + h\a(t k ,X k )\ + VhYT=i With, x k)€ik 
h\a(t k ,X k )\ + X^li "^,^)^ 



1 + h\a(t k ,X k )\ + V^IXi Wr(tk, X k )€ik\ 

m 

Kh\a(t k ,X k )\ + KVhJ2Wr{t k ,X k )\. 

r=l 

Using (USD and ([331), we obtain from (j3"H]) - (j3TT2"j) : 



< 



A < coh + c[\X k \ 2 h + KhJ2Wr(tk,X k )\' 



r=l 



h\a(t k ,X k )\ +Vhy)\a r (tk,Xk)\ 



r=l 



(3.13) 



+Kh 2 + K\X k \ 2 *h 2 
< Kh(l + \X k \ 2 + \X k \ 2 "h + \X k \^h 1 ' 2 ) < Kh(l + |X fc | 2 + \X k \^h 1 / 2 ). 



Now consider the last term in (13. 6h : 

E X ^ B[ .H|x fc | 2p - / |x fc+1 -x fc | i 



(3.14) 



< KE X ^ h (u)\X k \ 



2p-l 



h l \a(t k ,X k )\ l + h l / 2 J2Wr(U,X k )\% k \ 



r=l 



< KExQ R k (u) \X k \ 2p ~ l h l l 2 \l + \X k \ l 

where we used (|2.6p and (j3.3[) again as well as the fact that Xq r fe ( w ) an( ^ -^fc are ^t k -measurable 
while £ rfc are independent of Ft k ■ 

Combining (jg^j) . (IBTTjl . (IBTTUjl . and (JXHD, we obtain 



E Xn K , fc+1 (^)l^+i| 2p 
< E XClRk (u)\X k \ 2 P + JJ7iEx fljl >) I^! 2P " 2 [l + l^l 2 + \X k \ 3 ^/ 2 



(3.15) 



2/* 



+K^2E X ^ lRk (u J )\X k \ 2 ^ l h 1 / 2 1 + \X 



1=3 



2p 



< 



E XhR k {u)\X k \ 2 ? + KhExn (u) \X k \ 2p + Kj2®Xn Rk (") \X k \ 2p ~ l 



1=2 



+Kh 3 / 2 E XClRk (u;)\X k 



2p-2+3x 



2p 



+ Kh^Xn R . k ^)\X k \ 2p+l ^- 1) h l / 2 -\ 

1=3 



Choosing 



R = R(h)=h~ 1 /^- A \ 



(3.16) 



9 



we get E X ^» \X k \^^h^ < %w ( W ) \X k \^ and Xfl^^H l**!**^ 1 ' ^ < 
Xn Rh fc ( a; ) l^fc| 2p ; I = 3, . . . ,2p, and hence we re- write (|3. 15|) as 

* EXfl^HI^I* + KhE x ^ h Ju) \X k \^ + K±E Xflm (u>) \X k \^ h l 

1=1 

* E *n R(h J")\X k \ 2p + KMxa^Ju) + Kh, 

where in the last line we have used Young's inequality. From here, we get by Gronwall's inequality 
that 

E XClRwh (cu)\X k \ 2 P < K(l + E\X Q \ 2 P), (3.18) 

where R(h) is from (|3. 16[) and K does not depend on k and h but it depends on p. 
It remains to estimate ^X\ Rh k ( UJ )\^k\ 2p ■ We have 

X\ Rtk = 1 - Xn R k = 1 - Xn HjJk _ 1 X|jr Jk |<B = Xl^^ + Xn^^x^R 
k 

= ■■■ = y^Xn n ,_,X\x,]>r^ 

1=0 

where we put Xq r 1 = 1- Then, using (13. 4p . (I3.18p . (I2.1|) . and Cauchy-Bunyakovsky's and 
Markov's inequalities, we obtain 



1=0 

<(E|X + fe|^) 1/2 ^(E[ XWi|Xi|>i?(h) 

1=0 

=( e i^+^) i/2 eKwj^i>^: 



z=o 



* (r(y- |X| 2 ( 2 p +1 )( 6,< - 4 )^ 
< (E\X + k\^) 1/2 J2±— 



1=0 



i2 ( /i )(2 P +l)(6x-4) 



< K (E|X + k\ 4 P) l/2 (E(l + |Xo| 2 ( 2 p +1 )( 6 - 4 ))) 1/2 fc/i 2 ^ 1 < K(l + E |x | 4 f+ 2 ( 2 P+ 1 )( 6 ^- 4 )) 1 /2 ) 

which together with (|3.18p implies (|3.2p for integer p > 1. Then, by Jensen's inequality, 
holds for non-integer p as well. □ 

The next lemma gives estimates for the one-step error of the balanced scheme (|3.ip . 



Lemma 3.2. Assume that (|2.5p holds. Assume that the coefficients a(t,x) and o~ r (t,x) have 
continuous first- order partial derivatives int and that these derivatives and the coefficients satisfy 
inequalities of the form (|2.3p . Then the scheme (|3.ip satisfies the inequalities (|2.9p and (|2.10p 
with q± = 3/2 and q2 = 1, respectively. 
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The proof of this lemma is a routine analysis of the one-step approximation corresponding 
to (|3.ip using the equalities (|3.8|) - (|3.9|) . Since such analysis is similar to those done in the 
global Lipschitz case |13|, 1 15 j . we omit these routine calculations here. Lemmas 13.11 and 13,21 and 
Theorem 12.11 imply the following result. 

Proposition 3.1. Under the assumptions of Lemmas \3.1\ and \3.2\ the balanced scheme (|3.ip has 
mean-square order 1/2, i.e., for it the inequality (|2.13|) holds with q = q2 — 1/2 = 1/2. 

Remark 3.1. In the additive noise case the mean-square order of the balanced scheme (|3.ip 
does not improve (q\ and q2 remain 3/2 and 1, respectively). 

4 Fully implicit schemes 

Fully implicit (i.e., implicit both in drift and diffusion coefficients) mean-square schemes were 
proposed in p3] (see also [151 Chapter 1]), where their convergence was proved under global 
Lipschitz conditions. Here we analyze these schemes under the following assumptions, which 
are stronger with respect to the diffusion coefficient than Assumption 12.11 used in the previous 
Sections [2] and [3l 

Assumption 4.1. (i) The initial condition is such that 

E\X \ 2p <K<oo, for all p>l. (4.1) 

(ii) There exists a constant c\ > such that 

{x-y,a{t,x) -a{t,y)) < Cl \x-y\ 2 , t G [t , T], x, y G R d . (4.2) 

(Hi) There exist C2 > and yc > 1 such that 

\a(t,x)-a{t,y)\ 2 <c 2 (l + \x\ 2 "- 2 + \y\ 2 "- 2 )\x-y\ 2 , te[t ,T], x,yeR d . (4.3) 

(iv) The coefficients a r (t,x) have continuous bounded first-order spatial derivatives so that there 
are constants L\ > and L 2 > : 

\Va r (t,x)\<L 1 ,r=l,...,m, t G [t ,T], x G R d , (4.4) 

and 

\Va r (t,x)a r (t,x)-Va r (t,y)a r (t,y)\ < L 2 \x-y\, ,r = l,...,m, t G [t , T], x, y G R d . (4.5) 

In proofs which follow we will need some implications of Assumption 14.11 The condition 
(|4.2p implies that there is c > 

(z,a(t,x)) < c(l + |x| 2 ), t G [t ,T], x G R d . (4.6) 
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It follows from (fO|) that 

\a r (t,x)-a r (t,y)\<Li\x-y\, t G [t Q , T], x, y G R d , (4.7) 

and hence 

\a r (t,x)\ < Lx\x\ +L , (4.8) 
where Lq = max t6 r t0) rj \a r (t, 0)|. Further, there is L > : 

\Va r (t,x)a r (t,x)\ < L(l + \x\), t G [t , T], x G R d , (4.9) 

and 

|cr r (t,x)| 2 < L(l + |x| 2 ), i G [t , T], x £ R d . (4.10) 

For definiteness, we consider the following one-parametric family of methods for (jl.ip from 
the broader class of fully implicit schemes of [14} [15] : 

X k+l = X k + a(t k+x ,{l-X)X k + XX k+l )h (4.11) 

m d n 

- A EE£j(^' ( X - A ) Xfc + AX fc+1 )^ r (4 +A! (1 - \)X k + \X k+1 )h 

r=l j=l 



+ a r (t fc+A , (1 - X)X k + \X k+1 ) (C rh ) k Vh, 



r=l 



where < A < 1, t k+ \ = t k + Xh and (Cr/Jfc are i-i-d. random variables so that 

C h = * 

with £ ~ 7V(0, 1) and A h = y/2l\ \nh\ with I > 1. We recall [HI Lemma 1.3.4] that 



A, £ > A h , (4.12) 
A, £ < -A h , 



E(C 2 - (1) = (1 + 2VW^h\)h l - (4.13) 

Remark 4.1. Three choices of X are most notable: A = gives the explicit Euler scheme which 
is divergent JH [?|/ in the considered setting; A = 1 gives the fully implicit Euler scheme; and 
A = 1/2 corresponds to the mid-point rule, which in application to a system of Stratonovich SDE 
is derivative free jll5\, p. ^5]. 

Now we will study properties of the method (|4.1ip . 
Consider the one-step approximations corresponding to (]4.1ip 

m d n 

X = X x = x + a(t + Xh,U x )h- Xy2y2j^{t + Xh,U x )ai(t + Xh,U x )h (4.14) 

OX*' 

r=l j=l 

771 

+ Y,°r{t + Xh,U X )( rh Vh, 
r=l 
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where 



Note that 



U = U X := (l-X)x + XX > 



(4.15) 



m d 



U x = x 



+ Xa(t + Xh, U x ) - >?^2^2^f L (t + Xh, U x )4(t + Xh, U x )h (4.16) 



r=l i=l 

+A ^ a r (t + A/i, U x )C rh Vh. 



r=l 



Lemma 4.1. Lei < A < 1. Assume that Assumption ^. 1\ holds. For an arbitrary < e < 1, 
/md /iq > such that 



X 



hQCi + mXLih§ + mLi y2l^o| l n ^ol 



1 - s. 



(4.17) 



Then the equation (14. 14ft /or any < h < h$ has the unique solution X which satisfies the 
inequalities for some K > : 



\X-x\< K(l + \x\")h + K(l + |x|)V^| 



and 
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l^l 2 < ^(^o + l)VWM + ^[(l - A) 2 + ^M 2 , * e [*o, 21, s e 



Proof. Let 

a(t,x) = a(t,x) - X^2^2-^(t,x)al(t,x). 
r=l j=l 

For any fixed A, i, £ rft , and h, we introduce the function 



(4.18) 



(4.19) 



(4.20) 



i/j(z) = z — Xa(t + Xh, z)h — X a r (t + Xh, z)( rh Vh 

which is continuous in z due to our assumptions. The equation (|4.14p can be written as 

^(U x ) = x. (4.21) 

Using flJ2J), and (|477|) . we obtain 



(z - y,ip(z) - ip(y)) > \z - y\ 2 - hXci\z - y\ 2 - hmX 2 L 2 \z - y\ 



(4.22) 



-mXLi^/2lh\ \ah\\z - y\ 2 



= (1 - A /ici + mXL 2 h + mL x ^JHh\ \nh\ )\z - y\ 2 > e\z - y\ 2 > 0, 

i.e., ip(z) is uniformly monotone function for h < ho. This implies (see, e.g. |20|. Theorem 6.4.4, 
p. 167]) that (|4.14p has a unique solution. 
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We obtain from (ET2TD and (j4T22l) : 

e\U\ 2 < (U,il>(U)-MO)) = (U,x-i/,(0)) 

e,.„ 2 2, l2 2.,,„.., e lTTl2 2, ., 2A(L + l)v® In h\ 
< -\U\ 2 + -\x\ 2 + -|-0(O)| 2 < -\U\ 2 + -\x\ 2 + — ^ ^ 1 L, 

from which (|4.19p follows. 

|f/| 2 < Ja( A ( L o + 1)\/2^| In fcj + |x| 2 ) (4.23) 

which implies (|4.19|) . 

Further, it follows from (|4"35]) . P~2TJ) and (f¥722|) that 

m 

Ae|X-x| 2 = g|C7 - x| 2 < (U - x, -Xa(t + Xh,x)h - X^a r (t + \h,x)C, rh \fh) 

r=l 

< X 2 \X - x\ yi\a(t + Xh, x)\ + ^2^1 In h\ ^ + A/i, x)|J . 

Then, using (|4.3p and (|4.8p . we obtain (|4.18p . which completes the proof of Lemma 14. II for the 
implicit method (|4.1ip . □ 

Now we consider boundedness of moments of (|4.1ip . 

Lemma 4.2. Let 1/2 < A < 1. Assume that Assumption \4-l\ holds. Then for all < h < ho 
with ho from ()4.17p and for all k = 0, . . . , N the following inequality holds for the fully implicit 
scheme (|4.1ip for p > 1 : 

E|X fe | 2p < K(l + E|X | 2p ), (4.24) 

where K > is a constant. 

Proof. We note that f)4.6[) and fj4.9|) imply 

(x,a(t,x)) < (c + 3mAL/2)(l + |x| 2 ), i6[t 0) T], sGR, (4.25) 
which together with (|4.4p ensures that the solution of (jl.ip has all moments (|2.5p . p> 1 ^j. 
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Let £4+i = (1 - X)X k + \X k+1 (cf. (gJSD ). We have 

14 + i := |X fc+1 | 2 - |X fc | 2 = 2(U k+1 ,X k+1 - X k ) - (2A - l)\X k+1 - X k \ 2 (4.26) 

m 

= 2\h(U k+ i,a(t k+x ,U k+1 )) + 2\Vh(U k+1 ,^2a r (t k+x ,U k+1 ) (( rh ) k ) 

r=l 
m 

-(2A - l)/i 2 |a(i fc+A , ^ +1 )| 2 - (2A - l)h\ °r(t k+ x, U k+l ) (( rh ) k | 2 

r=l 

m 

-2(2A - l)h 3 / 2 (a(t k+x , U k+1 ), ^ <Tr(tk+\, U k+l ) (Crh)k) 

r=l 

m 

= 2X(U k+ i, ha(t k+x ,U k+1 )) + 2XVh(X k ,^a r (t k+ x, U k +i) (C r h)k) 

r=l 

m 

+(2A 2 - 2A + l)h\ J2<rr(t k+X , U k+l ) (( rh ) k | 2 - (2A - l)h 2 \a(t k+x , U k+1 )\ 2 

r=l 

m 

+2(1 - X) 2 h 3 / 2 {a(t k+x , U k+1 ), °r(t k+ \, U k+1 ) (C rh ) k ). 

r=l 

Expanding a r (t k+x ,U k+ i) at (t k+x ,X k ), we obtain 

m 

V k+1 = 2X(U k+1 ,ha(t k+x ,U k+1 )) + 2XVh(X k ,Y,*r(tk + \,X k ) (( rh ) k ) (4.27) 

r=l 

m 

+2XVh(X k ,J2V<Tr(tk+\,0)(U k+ i - X k ) (( rh ) k ) 

r=l 

m 

+(2A 2 - 2A + l)h\ ]T <rr(tk+x, U k+l ) (( rh ) k \ 2 - (2A - l)h 2 \~a(t k+x , U k+l )\ 2 



r=l 



+2(A - l) 2 h 3 / 2 (a(t k+ x, U k+1 ), Y °r(t k+ x, U k+1 ) (Crh) k ) 

r=l 



2X(U k+1 , ha(t k+x , U k+1 )) + 2XVh(X k ,J2 °r(t k+x , X k ) (C rh ) k ) 

r=l 

m 

+2X 2 h 3 / 2 {X k , Y V<7 r (i fc+A , 6)a{t k+x , U k+1 ) (( rh ) k ) 



r=l 



+2X 2 h(X k ,Y V^ r (t fc+A , 9) ]T <ri(t k+ x,X k ) (C lh ) k (C r fc) fc ) 

r=l 1=1 



rn 



+(2A 2 - 2A + l)h\ Y °r(tk+x, U k+1 ) (C rh ) k | 2 - (2A - l)h 2 \a(t k+Xl U k+1 )\ 2 

r=l 

m 

+2(1 - X) 2 h 3 / 2 {a(t k+x , U k+1 ), Y °r(t k+ x, U k+1 ) (C rh ) k ), 

r=l 

where 6 = vU k+ \ — (1 — v)X k , v G [0, 1], is an intermediate point. Using (|4.25p . (|4.10p . Young's 
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inequality, (|4.4p and (|4.19p . we obtain 

m 

V k+1 < Xh(2c + 3\mL)(l + \U k+1 \ 2 ) + 2\Vh{X k ,Y,Vr(tk+\,X k ){( rh ) k ) (4.28) 



r=l 



2A — 1 2A 4 m 

.^—h 2 \a(t k+x , U k+1 )\ 2 + — -h\X k \ 2 m ]T \Va r (t k+x , 9)\ 2 \ (( rh ) h 

r=l 



+\ 2 hm\X k \ 2 J2 Na r (t k+X , 0)\ 2 \ (C rh ) k | 2 + X 2 hmJ2 kl(*k+A, ^fc)| 2 | (Cifc)* l' 

r=l i=l 
m 

+ (2A 2 - 2A + l)fcm ^ M^ +A , U k+1 )\ 2 \ (( rh ) k \ 2 - (2A - l)h 2 \a(t k+x , U k+1 ) 

v) 4 m 

_ — hm J2Wr{t k+ X,U k+l )\ 2 \ (C rh ) k | 2 



r=l 



H ^ — /i 2 |a(i fc+A , t/ fc+ i)| 2 + 



r=l 



< A/i(2c + 3AmL)(l + \U k+1 \ 2 ) + 2\Vh(X k , J2 °r{tk+\, X k ) (C 



rh)k) 



r=l 



2A 



+ A2 [^3l + l]^?l^fcl 2 "»E I (Crfc)* I' + A 2 fcmL(l + \X k \ 2 ) £ | (C 



lh)k 



r=l 



Z=l 



+ [2A 2 _ 2A + 1 + ]M1 + l^+il 2 ) E I (C 



z/Jfc 



+ 2Av^(X fe ,^<7 r (t fc+A ,X fc )(C 



rh) k )i 



Then using (|4.23p . we arrive at 

(m 
l + El(Crh)* 
r=l / r=l 

where -KT > is independent of h and /c while it depends on A and on constants appearing in 
Thus 



l + |X fc+1 | 2 < l + |X fe | 2 +i^(l+|X fe | 2 ) 1 + E I (Crfc)fc I' ) +2AV^(X fe ,E^r(^+A,^) (Crfc)fc) 
Then for integer p > 1 we get 



r=l 



r=l 



[i + \x k+1 \ 2 y < (i + \x k \ 2 Y + K(i + \x k \ 2 Y52h i 



i=i 



1 + E I (Crfc)fc 



r=l 



+kE(i + I^IT"^ V2 



{=1 



(■Xfc, E °>(*fc+A, ^fc) (Cr/i)j( 



r=l 



whence, observing that X k are .7-^ -measurable while £ rk are independent of Tt k , it is not difficult 
to obtain 



E(l + |X fe+1 | 2 ) p < E(l + |X fc | 2 f + ™(l + |X fc | 



2\P 



\2\P 



+Kj2m(i + \x k \ 2 y- l h 1 / 2 

1=2 

< E(l + \X k \ 2 ) p + KhE(l + \X k \ 2 ) p , 
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(X k ,y^ j a r (t k+Xl X k ) (Cr/Jfc) 



r=l 



which together with Gronwall's inequality completes the proof of the lemma for integer p > 1. 
Then by Jensen's inequality for non-integer p > 1 as well. □ 

We have not succeeded in proving boundedness of moments for the mid-point scheme, i.e., 
(|4,lip with A = 1/2 under Assumption 14.11 One can observe that the proof of Lemma 14.21 is 
not applicable to this choice of A as the estimate in f|4.28j) blows up when A — > 1/2 and it is 
clear that the mid-point scheme is the boundary case. We also know [3] that for a r = ()4. is 
B-stable for A > 1/2 and not B-stable (in fact, not A-stable) for A < 1/2. It is natural to expect 
that for A < 1/2 the moments of (|4.1ip are not bounded and hence the method with A < 1/2 
is divergent under Assumption 14.11 (see also such a conclusion for the drift-implicit 8- method in 
|21|). In our experiments (Section [5]) the mid-point method produced accurate results. 

At the same time, we proved boundedness of moments for the mid-point scheme if in addition 
to Assumption 14.11 we require that the diffusion coefficients a r (t,x) are bounded. The proof is 
similar to the proof of Lemma 13.11 

Lemma 4.3. Let the assumptions of Lemma \4.2\ hold and in addition assume that the diffusion 
coefficients a r (t,x) are uniformly bounded. Then the moments of the mid-point method ([4. lip 
with A = 1/2 has bounded moments: for p > 1: 

E\X k \ 2p < K(l + E|X | 4(p+1)x - 4 ) 1/2 , (4.29) 

where K > is a constant. 

Proof. For k = 1 (cf. (|4.3|) ) . i.e., the global Lipschitz case, boundedness of moments of X k is 
established in (THUS]. Let h> 1. 
From (|4.26p . we have 

m 

\X k+ i\ 2 - \X k \ 2 = h(U k+ i,a(t k+1/2 ,U k+1 )) + \ /f h(U k+ i,J2 cr r(t k+1 /2,Uk+i) (Crh)k)- 

r=l 

Then using (|4.25p and Young's inequality, boundedness of a r and (|4.23p . we get 

m 

\X k+1 \ 2 - \X k \ 2 < h(c + 3\mL/2)(l + \U k+1 \ 2 ) + h\U k+1 \ 2 + KY,\(Crh) k \ 2 

r=l 

m 

< Kh(l + \X k \ 2 ) + KY,\(C rh ) k \ 2 , 

r=l 

from which one can obtain for integer p > 1 : 

E (1 + \X k \ 2 ) p < KhT p E (1 + \X \ 2 ) P . (4.30) 
Further, using (|4.25p . Young's inequality, boundedness of <7 r , (|4.4p . (|4.3|) and (|4.23p . we get from 
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m 

\X k+ i\ 2 - \X k \ 2 = (U k+1 ,ha(t k+ i /2 ,U k+1 )) + Vh(X k ,J2 a r(.t k+1/2 ,X k ) (( rh ) k ) 

r=l 

+-/ l 3 / 2 (x fc , v<7 r (t fc+ i, 0)a(t fc+1/2 , tr fc+1 ) (c rh ) fc ) 

Z r=l 

^ m m m 

+ VfJ r .(t fc+1 /2,6') ^V;(£ fc+ i/ 2 ,Xfc) (C^) fc (CrlOfc) + o^' ^ a r( t k+l/2, U k+1 ) (Cr/Jfc 

r=l Z=l r=l 

^ m 

+ -/i 3 / 2 (a(t fc+A , E/fc+l), ^2 a r(t k+ \, Uk+l) (C r , 



r=l 



r=l r=l 
5 rh)fc 

r=l 



< + |x fc | 2 )(i + ^ | (c rfc ) fc | 2 ) + Vh(x k , £ *r(t k+1/2 ,X k ) (c rh ) k ) 

r=l 

m 

+Kh^(l + \X k \^)Y,\(Cr 
r=l 

Choosing R(h) = /i" 1 / 2 ^" 1 ), we get after some additional calculation (cf. (|3.18p ): 

Exn R{h) k (u;)\X k \ 2 P < K(l + E\X \ 2 'P), (4.31) 

where £lR(h),k is the event as in (|3.5|) . 

Now using ()4.30p . (I4.3ip . and Cauchy-Bunyakovsky's and Markov's inequalities, we arrive at 

(cf. 

E XAjs{h) fc ( W )|X fc | 2 ^ < K(l + E|Xo| 4 ^ +1 )- 4 ) 1 / 2 , 

from which together with (|4.3ip the inequality (|4.29p follows. □ 

The next lemma gives estimates for the one-step error of the method (|4.11|) . 

Lemma 4.4. Let < A < 1. Assume that (|2.5p holds. Assume that the coefficient a(t,x) 
has continuous first order partial derivative in t and in x l and that the derivatives and the 
coefficient satisfy inequalities of the form (|4.3p ; the functions o~ r (t,x) have continuous first- 
order partial derivatives in t and that the derivatives and the coefficients satisfy inequalities of 
the form (|4.4p — (|4.5p ; and the functions V a r (t,x)o~ r (t,x) have continuous first partial derivatives 
in t and in x l which satisfy inequalities of the form (|4.5p . Then the method (|4.1ip satisfies the 
inequalities (|2.9p and (|2.10p with q\ = 2 and q% = 1, respectively. 

Proofs of this lemma is rather routine and similar to the global Lipschitz case |14l [15] and 
it is omitted here. Using Lemmas 14.1114.41 t ne next proposition follows from Theorem 12.11 

Proposition 4.1. Let for 1/2 < A < 1 the assumptions of Lemmas and \4-4\ hold and for 
A = 1/2 in addition assume that the diffusion coefficients a r (t,x) are uniformly bounded. Then 
the fully implicit method (|4.1ip has mean-square order 1/2, i.e., for it the inequality (|2.13p holds 
with q = 1/2. 
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Remark 4.2. Consider the commutative case, i.e., when Ajov = A r <7j (here the operator A r := 
(o~ r ,d/dx)) or in the case of a system with one noise (i.e., m = 1). Then in the setting of 
Lemma \4--4\ the mid-point method, i.e., (|4.1ip with A = 1/2, satisfies the inequalities (|2.9p and 
(I2,10p with qi = 2 and q% = 3/2, respectively (see such a result in the global Lipschitz case in 
\14\ \15Jl). Therefore, it converges in this case with mean-square order 1 when its moments are 
bounded. 



5 Numerical examples 

In this section we will test the following schemes: the balanced method (|3.ip from Section (3j the 
drift-implicit scheme (|2.17p ; the fully implicit Euler scheme (|4.1ip with A = 1; the mid-point 
method (|4.1ip with A = 1/2; the drift-tamed Euler scheme (a modified balanced method) [8]: 

(X \ m 

X k+l = X k + h a \ k + o- r (t k ,X k )Z rk Vh; (5.1) 
the fully-tamed scheme [TJ: 

a(X k )h + Y^T=i o- r (tk,X k )£ rk Vh 



Xk+l — X k H 

max ( 1, h 

and the trapezoidal scheme: 



ha{X k ) + YJ?=i °~r(t k , X k )^ rk y/h 



(5.2) 



h m 
X k+l = X k + - [a{X k+1 ) + a(X k )] + ^ <^r{t k , X k )£ rk Vh. (5.3) 

1 r=l 
As before, £ rk = (w r (t k+ i) — w r (t k ))/y/h are Gaussian Af(0, 1) i.i.d. random variables. We note 
that under Assumption 12.11 boundedness of second moments and strong convergence (without 
giving order) of ^-schemes, and in particular of (|5.3p . can be found in [2Tj. Strong convergence 
with order 1 /2 of (|5.ip under Assumption 14.11 is proved in [8] . Strong convergence of (|5.2p 
without order under Assumption 12 . 1 1 is proved in [7]. 

In all the experiments with fully implicit schemes, where the truncated random variables £ 
are used, we took I = 2 (see (|4.13p ). The experiments were performed using Matlab R2012a on a 
Macintosh desktop computer with Intel Xeon CPU E5462 (quad-core, 2.80 GHz). In simulations 
we used the Mersenne twister random generator with seed 100. Newton's method was used to 
solve the nonlinear algebraic equations at each step of the implicit schemes. 

We test the methods on the two model problems. The first one satisfies Assumption 14.11 

(nonglobal Lipschitz drift, global Lipschitz diffusion) and has two non-commutative noises. The 

second example satisfies Assumption 12. II (nonglobal Lipschitz both drift and diffusion). The aim 

of the tests is to compare performance of the methods: their accuracy (i.e., roughly speaking, 

size of prefactors at a power of h) and computational costs. We note that experiments cannot 

prove or disprove boundedness of moments of the schemes since experiments rely on a finite 

sample of trajectories run over a finite time interval while blow-up of moments in divergent 

methods (e.g., explicit Euler scheme) is, in general, a result of large deviations I16j. 
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Example 5.1. Our first test model is the Stratonovich SDE of the form: 

dX = (1 -X 5 )dt + Xo dwi + dw 2 , X(0) = 0. (5.4) 

In Ito's sense, the drift of the equation becomes a(t,x) = 1 — x 5 + x/2. Here we tested the 
balanced method (|3.ip ; the drift-tamed scheme (|5.ip ; the fully implicit Euler scheme (|4.1ip with 
A = 1; the mid-point method (|4.11|) with A = 1/2. We note that for all the methods tested on 
this example except the mid-point rule mean-square convergence with order 1 /2 is proved either 
in earlier papers or here as it was described before. 

To compute the mean-square error, we run M independent trajectories X^(t), x9^: 

(M \ 1/2 

l]T[X«(r)-X$] 2 j . (5.5) 

We took time T = 50 and M = 10 4 . The reference solution was computed by the mid-point 
method with small time step h = 10~ 4 . It was verified that using a different implicit scheme for 
simulating a reference solution does not affect the outcome of the tests. We chose the mid-point 
scheme as a reference since in all the experiments it produced the most accurate results. 

Table [5TT1 gives the mean-square errors and experimentally observed convergence rates for the 
corresponding methods. We checked that the number of trajectories M = 10 4 was sufficiently 
large for the statistical errors not to significantly hinder the mean-square errors (the Monte 
Carlo error computed with 95% confidence was at least 10 time smaller than the reported mean- 
square errors except values for (|5.ip at h = 0.1 and 0.05 where it was at least 5 time smaller 
than the mean-square errors). In addition to the data in the table, we evaluated errors for (|3.ip 
for smaller time steps: h = 0.002 - the error is 9.27e-02 (rate 0.41), 0.001 - 6.86e-02 (0.44). 
The observed rates of convergence of all the tested methods are close to the predicted 1/2. For 
a fixed time step h, the most accurate scheme is the mid-point one, the less accurate scheme 
is the new balanced method (|3.ip . To produce the result with accuracy ~ 0.06 — 0.07, in our 
experiment of running M = 10 4 trajectories the scheme (|5.ip required 170 sec, the mid-point 
(I4TTTD with A = 1/2 - 329 sec, (|4TTT1) with A = 1 - 723 sec, and (HEED - 1870 sec. That is, 
our experiments confirmed the conclusion of [8] that the drift-tamed (modified balance method) 
(|5.ip from [8] is highly competitive. We note that (|5.ip is not applicable when diffusion grows 
faster than a linear function and that in this case the balanced method (|3.ip can outcompete 
implicit schemes as it is shown in the next example 

Example 5.2. Consider the SDE in the Stratonovich sense: 

dX = (1 -X 5 )dt + X 2 o dw, X(0) = 0. (5.6) 

In Ito's sense, the drift of the equation becomes a(t, x) = 1 — x 5 + x 3 . 

Here we tested the balanced method (|3.ip ; the fully-tamed Euler scheme (|5.2[) ; the drift- 
implicit scheme (|2.17p ; the fully implicit Euler scheme (|4,lip with A = 1; the mid-point method 
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Table 5.1: Example \5.1[ Mean-square errors of the selected schemes. See further details in the 
text. 



h 


(14.111). A = 1 rate 


(I4.11j). A = 1/2 rate 


(pTT|) rate 


(I3TTD rate 


0.1 


1.712e-01 


1.443e-01 


3.748e-01 


3.594e-01 


0.05 


1.234e-01 0.47 


9.224e-02 0.65 


2.103e-01 0.83 


3.017e-01 0.25 


0.02 


7.692e-02 0.52 


5.261e-02 0.61 


9.472e-02 0.87 


2.297e-01 0.30 


0.01 


5.478e-02 0.49 


3.549e-02 0.57 


6.104e-02 0.63 


1.778e-01 0.37 


0.005 


3.935e-02 0.48 


2.487e-02 0.51 


3.959e-02 0.62 


1.354e-01 0.39 



(|4.1ip with A = 1/2; and the trapezoidal scheme f)5 . 3|) . We recall that in the case of nonglobal 
Lipschitz drift and diffusion, for the drift-implicit scheme (|2.17p and the balanced method f)3. 1 j) 
mean-square convergence with order 1/2 is shown earlier in this paper; strong convergence of 
the trapezoidal scheme (|5.3|) without order is proved in [21], it is natural to expect that its 
mean-square order is 1/2 which is indeed supported by the experiments. Strong convergence 
of (|5.2p without order is proved in [7]. We note that it can be proved directly that implicit 
algebraic equations arising from application of the mid-point and fully implicit Euler schemes 
to (|5.6p have unique solutions under a sufficiently small time step. 

The reference solution was computed by the mid-point method with small time step h = 10 -4 . 
The time T = 50 and M = 10 4 in (153]) . 

The fully-tamed scheme (|5.2p did not produce accurate results until the time step size is at 
least h = 0.005 and we do not then report its errors here but see the remark below. 

Remark 5.1. The fully-tamed scheme (|5.2|) appears to be of a low practical value. If at a step k*, 
the event O := ha(Xk) + Y^r=i a r{tk, Xk)£, r kVh > 1/h happens, then in the case of (|5.6p the 
trajectory X^, k > oscillates approximately between Xk, and Xk t — signum(Xk t )/h. Since 
the probability of the event O is positive for any step size h > and grows with integration time, 
it is unavoidable that in some scenarios (i.e., on some trajectories) such oscillatory behavior 
will appear. For instance, in this experiment for h = 0.1 we observed 989 out of 1000 paths for 
which O happened over the time interval [0,50]; for h = 0.05 - 866 out of 1000 paths. From 
the practical point of view, (|5.2p works as long as the explicit Euler scheme works ( cf. \10^ and 
also 1151 p. 17]). The strong convergence (without order) of (|5.2|) FT]/ in comparison with the 
explicit Euler scheme is due to the following fact. When event O happens for the Euler scheme 
its sequence X^ starts oscillating with growing amplitude which leads to unboundedness of its 
moments and, consequently, its divergence in the mean-square sense. For (|5.2p . the oscillations 
are bounded by ~ 1/h and since the probability of O over a finite time interval decreasing with 
decrease of h, then the moments are bounded uniformly in h. At the same time, the one-step 
approximation of (|5.2p does not satisfy the conditions (|2.9p and (|2.10p of Theorem \2.1\ We note 
that the explicit balanced-type scheme (|3.ip does not have such drawbacks as (|5.2D , 
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Tabic [531 gives the mean-square errors and experimentally observed convergence rates for the 
corresponding methods. We checked that the number of trajectories M = 10 4 was sufficiently 
large for the statistical errors not to significantly hinder the mean-square errors (the Monte Carlo 
error computed with 95% confidence was at least ten time smaller than the reported mean-square 
errors). In addition to the data in the table, we evaluated errors for (|3.ip for smaller time steps: 
h = 0.002 - the error is 3.70e-02 (rate 0.41), 0.001 - 2.73e-02 (0.44), 0.0005 - 2.00e-02 (0.45), i.e., 
for smaller h the observed convergence rate of (|3.1|) becomes closer to the theoretically predicted 
order 1/2. Since (|5.6p is with single noise, Remark 14.21 is valid here which explains why the mid- 
point scheme demonstrates the first order of convergence. The other implicit schemes show the 
order 1/2 as expected. Table [531 presents the time costs in seconds. Let us fix tolerance level at 
0.05 — 0.06. We highlighted in bold the corresponding values in both tables. We see that in this 
example the mid-point scheme is the most efficient which is due to its first order convergence in 
the commutative case. Among methods of order 1/2, the balanced method (|3.ip is the fastest 
and one can expect that for multi-dimensional SDE the explicit scheme (|3.ip can considerably 
outperform implicit methods (see a similar outcome for the drift-tamed method (|5.ip supported 
by experiments in [8]; note that (|5.ip . in comparison with (|3.ip . is, as a rule, divergent when 
diffusion is growing faster than a linear function on infinity). 



Table 5.2: Example \5.2l Mean-square errors of the selected schemes. See further details in the 
text. 



h 


(|2.17p rate 


(I4.1ip. A = 1 rate 


<!4.1ip. A = 1/2 rate 


rate 


(fBTTj) rate 


0.2 


3.449e-01 


1.816e-01 


1.378e-01 


4.920e-01 


2.102e-01 


0.1 


2.441e-01 0.50 


1.331e-01 0.45 


8.723e-02 0.66 


3.526e-01 0.48 


1.637e-01 0.36 


0.05 


1.592e-01 0.62 


9.619e-02 0.47 


5.344e-02 0.71 


2.230e-01 0.66 


1.270e-01 0.37 


0.02 


8.360e-02 0.70 


6.599e-02 0.41 


2.242e-02 0.95 


1.048c-01 0.82 


9.170e-02 0.36 


0.01 


5.460e-02 0.61 


4.919e-02 0.42 


1.145e-02 0.97 


5.990e-02 0.81 


7.065e-02 0.38 


0.005 


3.682e-02 0.57 


3.522e-02 0.48 


5.945e-03 0.95 


3.784e-02 0.66 


5.393e-02 0.39 



Table 5.3: Example \5.S\ Comparison of computational times for the selected schemes. See 
further details in the text. 



h 


(I2.17p 


(14.111). A = 1 


(|4.1ip. A = 1/2 






0.2 


9.25e+00 


1.10e+01 


9.33e+00 


1.20e+01 


3.98e+00 


0.1 


1.77e+01 


2.17e+01 


1.80e+01 


2.30e+01 


7.49e+00 


0.05 


3.42e+01 


4.26e+01 


3.51e+01 


4.48e+01 


1.41e+01 


0.02 


8.33e+01 


1.04e+02 


8.69e+01 


1.10e+02 


3.37e+01 


0.01 


1.64e+02 


2.05e+02 


1.73e+02 


2.19e+02 


6.62e+01 


0.005 


3.25e+02 


4.07e+02 


3.47e+02 


4.37e+02 


1.32e+02 
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A Proof of the fundamental theorem 

Note that in this and the next section we shall use the letter K to denote various constants 
which are independent of h and k. The proof exploits the idea of the prove of this theorem in 
the global Lipschitz case |llj . 

Consider the error of the method Xt ,x (tk+i) at the (k + l)-step: 

p fe+ i := X t0)Xo {tk+i) ~ X t0tXo (tk+i) = X tk)X (t k )(tk+i) - X tkt x k (tk+i) (A.l) 
= i x t k ,x(t k ){tk+i) ~ X tk:Xk (tk+i)) + (X tk ,x k (tk+i) ~ X tki x k (tk+i)) ■ 
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The first difference in the right-hand side of (jA.ip is the error of the solution arising due to the 
error in the initial data at time t k , accumulated at the A:-th step, which we can re- write as 

s t k ,x(t k ),x k (tk+i) = S k+ i := X tk:X (t k )(tk+i) - X tk)Xk {tk+i) = Pk + z t k ,x(t k ),x k {tk+i) 
= Pk + z k+l, 

where Z is as in (|2.14p . The second difference in (jA.ip is the one-step error at the (k + l)-step 
and we denote it as r k+ \ ■ 

rk+i = X tk)Xk {tk+i) ~ X tkt x h (t k +i). 
Let p > 1 be an integer. We have 

E\p k+1 \ 2p = E|5 fc+1 + r fc+1 | 2p = E[(5 fc+1 ,5 fc+1 ) + 2(5 fc+1 ,r fc+1 ) + (r fc+1 ,r fc+1 )]P (A.2) 

< E|S fc+1 | 2p + 2pE | S k+1 \ 2p - 2 (p k + Z k+1 ,r k+1 ) + KY,^\Sk + i\ 2p ~ l \rk+i\ l . 

1=2 

Due to (|2.15p of Lemma I2.1| the first term on the right-hand side of (|A.2|) is estimated as 

E\S k+1 \ 2p <E\p k \ 2p (l + Kh). (A.3) 

Consider the second term on the right-hand side of (|A.2|) : 

E|S, +1 | 2p " 2 (p k + Z k+1 ,r k+1 ) = E\p k \ 2p ~ 2 ( Pk ,r k+1 ) (A.4) 
+E (\S k+1 \ 2p ~ 2 - \ Pk \ 2p ~ 2 ) (p k ,r k+1 ) + E\S k+1 \ 2p ~ 2 (Z k+U r k+1 ). 

Due to J-i fe -measur ability of p k and due to the conditional variant of (|2.9p . we get for the first 
term on the right-hand side of (|A.4p : 

E \ Pk \ 2p ~ 2 ( Pk ,r k+1 ) < KE \p k \ 2p " 1 (1 + \X k \ 2a ) 1 ' 2 h^. (A.5) 

Consider the second term on the right-hand side of (|A.4p and first of all note that it is equal to 
zero for p = 1. We have for integer p > 2 : 

2p-3 

E (\S k+1 \ 2p ~ 2 - \ Pk \ 2p ~ 2 ) ( Pk ,r k+1 ) < KE\Z k+1 \ \ Pk \\r k+1 \ £ \S k+1 \ 2p ~ 3 ~ l \p k \ l . 

1=0 

Further, using J-i fe -measur ability of p k and the conditional variants of (|2.10p . (|2,15p and (|2.16f) 
and the Cauchy-Bunyakovsky inequality (twice) , we get for p > 2 : 

E(\S k+1 \ 2p ~ 2 -\p k \ 2p - 2 )(p k ,r k+l ) (A.6) 

< KE |p fc | 2p_1 (1 + \X{t k )\ 2H ~ 2 + |X fc | 2 — 2 ) 1 /4/ 1 92+i/2 (1 + \ Xk \ 2a ) l l 2 . 

Due to J 7 t fe -measurability of p k , the conditional variants of (|2.10p and (|2.16p and the Cauchy- 
Bunyakovsky inequality (twice), we obtain for the third term on the right-hand side of (|A.4p : 

E|S fc+1 | 2p - 2 (Z fc+1 ,r fe+1 )<E^ (A.7) 

< KE^" 1 (1 + \X(t k )\ 2 *- 2 + \X k \ 2K - 2 ) l l i h^ +1 ' 2 {l + \X k \ ia ) l '\ 
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Due to J-"i fc -measur ability of p k and due to the conditional variants of (|2.10p and (|2.15p and 
the Cauchy-Bunyakovsky inequality, we estimate the third term on the right-hand side of (|A.2|) : 

2p 2p 

Kj2E\S k+1 \ 2p - l \r k+1 \ l < ir^E[E(|5 fe+1 | 4 ^|^J 1 /2 E( | rfe+1 |2«| Ji ji/2 ] (A8) 

1=2 1=2 

2p 

< Kj2 E l\Pk\ 2p ~ l h lq2 {l + \X k \ 2la f' 2 }. 

1=2 

Substituting ([A.3j) - (|A.8jl in (|A.2[) and recalling that q\ > qi + 1/2, we obtain 
E|p fc+1 | 2p < E|p fc | 2p (l + Kh) + KE |p fc | 2p_1 (1 + \X k \ 2a ) l ' 2 h q2+l ' 2 

+KE\p k \ 2p ~ 1 (1 + \X(t k )\ 2 *- 2 + \X k \ 2x ~ 2)1/4^2+1/2^ + \x k \ 2a )V 2 

+KE\ Pk \ 2p ~ l (1 + \X(t k )\ 2 "~ 2 + |^ fc | 2 ^- 2 ) 1/4 ^ g2+1/2 (l + (Xfcl 40 ) 1 / 4 

2p 

+ir^E[| /3fe | 2p - / /i /92 (l + \X k \ 2al ) 1 / 2 } 

1=2 

< E\p k \ 2p {l + Kh) + KE IPkl 2 ^ 1 (1 + \X(t k )\ 2 *~ 2 + \X k \ 2x - 2 ) 1 / i h q2+1/2 {l + \X k \ 2a ) 1 ' 2 

2p 

+KY,n\Pk\ 2p ~ l h lq2 (i + \x k \ 2ia ) 1 / 2 ]. 

1=2 

Then using Young's inequality and the conditions (|2.5p and (|2.12|) . we obtain 

E|p fc+1 | 2p < E\p k \ 2p + KhE\p k \ 2p + K(l + E|X |^^- 1 ) +2pQ ^)/i 2p ( g2 - 1 /2)+i 

whence ()2. 13j) with integer p > 1 follows by application of Gronwall's inequality. Then by 
Jensen's inequality (|2.13p holds for non-integer p as well. □ 



B Proof of Lemma 12.11 

Lemma 12. II is an analogue of Lemma 1.1.3 in |15j . 

Proof. Introduce the process St xy (s) = S(s) := X tx (s) — X ty (s) and note that Z(s) = 
S(s) — (x — y). We first prove (j2. 15|) . Using the Ito formula and the condition (|2.2[) (recall that 
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2J2]) implies fl23J)), we obtain for 9 > : 



rt+e 

E\S(t + 9)\ 2p = \x- y\ 2p + 2p J E\S\ 2p ~ 2 



ST(a(t,X t , x (s))-a(t,X t , y (s))) 



+ *^Wr(t,XtM)-<rr(t,Xt, v (s))\' 



r=l 



ds 



rt+e 

+2p(p-l) y E|5| 2p " 4 
/•t+e 

< |x-y| 2 f + 2p y E|S| 2p " 2 

+-^^— / ^Mi,X^( S ))-a r (i,X^(s))| 
< \x -y\ 2p + 2p Cl j E\S{s)\ 2p ds 



ds 



S\a(t,X t , x ( S ))-a(t,X t)y (s))) 

ds 



from which (|2.15|) follows after applying Gronwall's inequality. 

Now we prove f|2. 16|) . Using the Ito formula and the condition (|2.2p . we obtain for 9 > : 



rt+e 

E\Z{t + 9)\ 2p = 2p J E\Z\ 2p ~ 2 



ZT(a(t,X t , x (s))-a(t,X t Js))) 



(B.l) 



^ m 

+ V^Wr(t,X t M) - a r (t,Xt,y(8))\' 



r=l 



ds 



+2p(p - 1) j' +<) E\Z\ 2p ~ A 



Zi J> r (i, X t , x (s)) - a r (t, X t>y (s))] 



r=l 



ds 



rt+e 

< 2p y E|Z| 2p " 2 (s) 



2b - i r t+e 171 

/ ^k r (t,x fiX ( s ))-a r (t,x^( s ))| 

2 ^ r=l 
/■t+0 

-2p y E|Z| 2 P- 2 (x-y,o(t,X tiX ( S ))-a(t,X^( S )))d S 



ST(a(t,X t , x (s))-a(t,X tiy (s))) 



ds 



rt+V rt+ti 

<2p Cl / E\Z\ 2p ~ 2 \S\ 2 ds-2p / E|Z| 2p - 2 (x-y,a(t,X M ( S ))-a(t,X^( S )))ds. 

Using Young's inequality, we get for the first term in the right-hand side of (|B,ip : 

rt+e r t+e 
2p Cl E\Z\ 2p ~ 2 \S\ 2 ds < 4pci / E\Z\ 2p ~ 2 (\Z\ 2 + \x-y\ 2 )ds (B.2) 
Jt Jt 

rt+e rt+e 

<K E\Z\ 2p ds + K\x-y\ 2 E\Z\ 2p ~ 2 ds. 
Jt Jt 
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Consider the second term in the right-hand side of (jB.ip . Using Hoelder's inequality (twice), 
flgHgp and (^3]) . we obtain 



-2pJ E\Z\ 2 P- 2 (x-y,a(t,X t , x (s))-a(t,X t!y (s)))ds (B.3) 
pt+e 

< 2pJ^ E\Z\ 2p - 2 \a(t,X t , x (s))-a(t,X tty (s))\\x-y\ds 

< K\x - y\ J* mz^] 1 - 1 ' 9 [E\a(t,X t , x (s)) - a(t, X t>y (s))\ p ] 1/p ds 



-t+e 



< K\x-y\J [E\Z\' 
x(E[(l + |X M .( S )| 2 - 2 + \X t , y (s)\ 2 «- 2 yl 2 \X t , x (s) - X tiy (s)\^ds 

< K\x - y\ |' + [E\Z\ 2p ] (E[(l + \X ttX (s)\ 2 *- 2 + \X t ^ S )\ 2 ^ 2 ) p }f ,2p 
x (E[\X t M - X t , y (s)\ 2p ]) 1/2p ds 

< K | a; _y|2 (1 + | ;r |2x- 2 + | y |2x- 2) l/2 ^ [E|Z|*] 1_1/P cfe. 

Substituting (|BT2|) and (fB~3l) in (fB~Tj) and applying Hoelder's inequality to E\Z\ 2p ~ 2 ■ 1, we 

get 

E|Z(i + 0)| 2p < K / E|Z| 2p ds + K|x-y| 2 (l + |x| 2>< - 2 + |y| 2x - 2 ) 1/2 / [El^l 2 ? 5 ] 1 1/P ds 

(B.4) 

whence we obtain ()2. 16j) for integer p > 1 using Gronwall's inequality as, e.g. in |19t p. 360], 
and then by Jensen's inequality for non-integer p > 1 as well. □ 
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